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INTRODUCTION 

This article gives a brief summary of some results obtained by Nasser (ref. 1) on modeling 
and simulation of inequality problems in multibody dynamics. In particular, the augmented 
Lagrangian method discussed here is applied to a constrained motion problem with 
impulsive inequality constraints. A fundamental characteristic of the multibody dynamics 
problem is the lack of global convexity of its Lagrangian. The problem is transformed into a 
convex analysis problem by localization (piecewise linearization), where the augmented 
Lagrangian has been successfully used [see Glowinski and Le Tallec (ref. 2); Glowinski, 
Lions, and Tremolieres (ref. 3); and Fortin and Glowinski (ref. 4)]. A model test problem is 
considered, and a set of numerical experiments is presented (Figures 1 through 9). 

MATHEMATICAL MODEL 


Functional Context 

X = H m (o, T; R") , (1) 


H m (0,T;R N } = {u: v t C m ~ 1 ([0, 77; R N ), € L 2 (o, T; R*) j , (2) 

£ is a Lagrangian function, J is a nonlinear functional, and K is a closed subset of X. 
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Main Problem 


Find 0 S K, for which J is stationary 


K = K 1 nK 2 , (3) 

K l = ^uiX:g.(u(t))=0 , j = 1,2 k a.e. on (0, Tjj , (4) 

K 2 = |t>S X:g.(v(t)) SO , j = k + l,...,k+l a.e. on (0, T) j . (5) 

The functions gftXt)) are real valued; finally, J is defined by 

C T 

J(v) = £(v,v,t)dt. 

j 0 


The stationarity of J at 0 can be formulated as shown in the following section. 

AUGMENTED L AGRAN GIAN FORMULATION 


Following a well-known technique [Glowinski and Le Tallec (ref. 2); Glowinski, Lions, 
and Tremolieres (ref. 3); and Fortin and Glowinski (ref. 4)], we associate to (l)-(6) the 
following problem: 

Find 0 and A, with QtK and AS A, for which the following augmented functional J r is 
stationary: 


J r (v,n) 


T 

£ (v, v, n) dt , 

0 


(7) 


where 


£ r (v,n) = £(v, v) + <n l ,G 1 (u)> + <n r GJv)> + 


2 [ <G r G l > + <G 2 ’ G 2 > } ’ 


( 8 ) 


G(v) = 


g/vl.gjv), 



(9) 
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There exists a large amount of literature dealing with the case K 2 = 0, which leads to 
index 3 differential algebraic equations. The case when K 2 is nonempty is considerably 
more difficult from a mathematical point of view, and hence fewer technical papers have been 
devoted to it. The methodology we shall present includes treatment of both cases. 

SOLUTION ALGORITHMS 

Given 9 k (tJ, B k (tJ, and X k (tj, compute B k+1 (tJ, B k+I (tj , and X k+I (tJ via the 
following: 

Vr(»*T®» + l- *»)=*’ < 14 > 

\*,=' , aK + o c (M , .))1 ; (15) 

P A is the projection operator associated with the set A. For the choice of r and p, see 
Nasser (ref. 1) and Glowinski, Lions, and Tremolieres (ref. 3). 

LINEARIZATION AND TIME DISCRETIZATION 

Following Nasser (ref. 1), we introduce the perturbation SB, SB, SB of B, 8, 6 to obtain the 
following system: 

M(B ) SB + C(B,B) SB + K(B,6 ) SB + R(B) SB -h S(B) SB = r.h.s. (B) . (16) 
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( 17 ) 


*+/ 

s = X^*/ w - 


1=1 


I k , , * + J / 

iK«/^*/9MvH 2 )+ X K«/<w ■*/»> + [ v,* w 

,=i v ' t=k+l 


(18) 


r.h.s.(d) = - 

k+l 

X 

l=k+l 


k+l 


XVW e) + X l,V> , + r X »/«•*»*/»> 


»=i 


l=k+l 


1=1 


k + l j 

+ x ■%*,*« 


(19) 


Taylor Series Expansion of 0 and 0 


Using Taylor series expansion for 0 and 0, we get 


0a + d<) = 6(0 + 6 (v (t)At + 6 <2) (0 + d (3) (0 + 6 <4> (0 + o(jf 5 ) , 


( 20 ) 


At 2 , M 


d(t + At) = e (1) (t + At) = e <u (t) + e (2) (t)At + e (3) (o — + e (4 \t) — + o(At 4 ) , (21) 

2 n ' * 


where dots or superscripts denote the order of the derivatives. 


Let 


86(0 = 6(t + At) - 6(t) , 


86(0 = 6(t + At) - 6(t) , 


(22) 

(23) 


6 (3) (0 


e <2) (t + At) - 6 <2) (0 


At 


(24) 
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66(0 = 6(t + At) - 6(0 , 


(25) 


se <2) = ef 3) (t)At. 


(26) 


Time Discretization of the Differential Equation (16) 

Linear Acceleration Method 

This is a widely used scheme in structural dynamics. It consists of assuming that terms 
involving d (4> (t) in equations (20) and (21) are negligible and that the acceleration between t 
and t + At varies linearly [i.e., according to equation (26)]. Substituting equation (26) into 
equation (20), we get 


S 6 

66 = — -66(0 - — 6(0 - 36(0 . (27) 

At 2 

Taking equations (26) and (27) into account in equation (21), we obtain 

56 = —5d(t) - 36(0 - —6(0 . (28) 

At 2 

Substituting equations (27) and (28) into equation (16), we get the following linear system 
(in 66): 


A66 = b , 


(29) 


where 


A 


6 3 

= —M+ —C-hK-hR +S, 
At 2 ^ 


(30) 


b = 


— 6(t)M -h 3M6(t) + 3C6(t) + —6(0C. 
At 2 


(31) 
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Higher Order Time Discretization Schemes 


We assume that terms involving derivatives of order 5 and higher are negligible and that 


e (3> (t) = 


e (2) (t + At) - e^a - ao 

2At 


( 2 ), 


(32) 


e t4) (o = 


e rz, (t + At) - 26 <2) (0 + & z/ (t - At) 
A? 


42 ), 


(33) 


and 

86(0 = 6(t + At) - 6(0 . (34) 

Substituting equations (32), (33), and (34) into (20)-(21) and rearranging terms, we get 
(analogous to the linear acceleration method): 

AS6 = b, (35) 


where 


8 10 

A = — M+ C + K + R+S, 

At 2 3At 


(36) 


b = 


-Mdj — Cd 2 , 


(37) 



LXl «• 

— 6(t - At) - 6(0 At - 
24 


13 At 2 
24 



(38) 


*2 = 


13At 

12 


.. At •• 5At 

e(t )--e(t-Ao+—d I 


(39) 


A second algorithm is as follows: Given 6(tJ, 6(tJ, 6(tJ, 6(t nI ), and A (tj, update 

“ d h+M via 
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(40) 


S6 k = A~*b , 


* k + l P A *^k +pG {^k+l)}’ 


(41) 


where 


A = 


R 


k+l 




R 


i — 1 1 2, . . . , k , 0 , i — k + l,.. .,k+l\ 


In fact. 


(42) 


A =L 2 (0,T;A*y 


(43) 


Algorithm (40)-(41) can be used if A and 6 are replaced by A* and b* . 

Other integration schemes, such as the ones in Dean, Glowinski, Kuo, and Nasser (ref. 5), 
may be used, also. 

The acceleration 6(0 may be updated from the solution of equation (14) after convergence 
on (S6(0,M0) has been achieved. 

Choice of r and p 

The parameters r, p, and At are the variables controlling the stability. For optimal 
choice of these parameters, refer to Glowinski and Le Tallec (ref. 2) and Glowinski, Lions, 
and Tremolidres (ref. 3). 

Using the projection method substantiated and systematically developed in Nasser 
(ref. 1), the equations of motion of the unconstrained system can be obtained in the following 
form: 

M8Q -hCSQ-h K8Q = 0 . (44) 

The projection method without piecewise linearization has been used by Keat (ref. 6) and is 
equivalent to the well-known Kane’s method (ref. 7). 
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TEST PROBLEM 


Consider a planar two-body system with a rigid obstacle, as shown in Figure 1. The 
Cartesian coordinates are related to the Lagrangian coordinates by 


x t = dj sin d { , 


(44) 

x 2 = lsinQ { + a 2 sind 2 , 


(45) 

y t = cosdj , 


(46) 

y 2 = lcosd 1 + a 2 cos & 2 , 


(47) 

2 KE = [m,(x, + y,) +/,«?] + [™ 2 (i 2 + j 2 ) 


(48) 


PE = rrijga^l - cosd^ + m 2 g — cosfl^ + a^l — cosd^j j 
The stationarity of the Lagrangian £ is given by 

(/,+ mfj + m 2 l)e i + m 2 la 2 cos(d 2 - - m^sin (fl, - 0,) + m J g/sm0 / 

+ m 2 glsind I = 0 , 


(49) 


(50) 


m 2 a $2 + m J a 2 COS ( 6 2 ~ d l)^l ~ m 2 la 2 


fi sin ( e 2- d l) 


+ m 2 ga 2 sind 2 = 0 . 


(51) 


Data: 

nij = mass of body 1 
m 2 = mass of body 2 
l = length 

I { = moment inertia of body 1 
I 2 = moment inertia of body 2 
g = acceleration of gravity 
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Constraints: 

gl = h s * n + d 2: 0 

%2 ~~~ & l ^ ^2 @2 ^ ^ 

For the case r = 0 , the augmented Lagrangian method reduces to the multiplier 
method used for the treatment of Coulomb or dry friction problems in Dean, Glowinski, Kuo, 
and Nasser (ref. 8). For the case r* 0 , 1=0, the scheme reduces to the well-known penalty 
method. The parameter r is the spring stiffness coefficient used in classical contact 
problems. 


CONCLUSION 

The augmented Lagrangian method successfully applies to contact/constrained motion 
problems of multibody dynamics. For constraints involving 0, the technique still applies; 
however, the details are rather lengthy and were omitted. The case of elastic bodies offers no 
mathematical difficulty except in the details, and the convergence is influenced by the 
spatial discretization largest mesh size. For further details, refer to Kikuchi and Oden 
(ref. 9) and Nasser (ref. 1). 
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Figure 6. Penalty method/high-order scheme - energy. 
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Figure 7. Penalty method/high-order scheme - constraint force 
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Figure 8. Augmented Lagrangian energy comparison for high-order scheme versus 
the linear acceleration method. 
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Figure 9. Penalty/explicit Euler scheme - energy. 
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